Time Series Analysis: Car Sales In Italy

# Load Library

library(data.table) # transpose
library(forecast) # base
library(ggplot2) # base
library(TSstudio) #ts_plot, ts_heatmap, ts_seasonal
library(tseries) # adf.test
library(tsutils) # BIC
library(funtimes) # notrend_test
library(flexmix) # mseastest

# Functions

my_plot.decomposed.ts = function(x, title="", ...) {
  xx <- x$x
  if (is.null(xx)) 
    xx <- with(x, if (type == "additive") 
      random + trend + seasonal
      else random * trend * seasonal)
  plot(cbind(observed = xx, trend = x$trend, seasonal = x$seasonal, random = x$random), 
       main=title, ...)
} # In order to add personal title

my_plot.lags.ts = function (ts.obj, lags = 1:12, main = "", margin = 0.02, Xshare = TRUE, 
    Yshare = TRUE, n_plots = 3) 
{
    `%>%` <- magrittr::`%>%`
    df <- df_wide <- p <- obj.name <- lag <- lag_plots <- time <- NULL
    obj.name <- base::deparse(base::substitute(ts.obj))
    if (!is.numeric(lags)) {
        warning("The 'lags' parameter is not valid, using the defualt setting (lags = 1:12)")
        lags <- 1:12
    }
    else if (base::any(lags <= 0)) {
        warning("The 'lags' parameter is not valid, using the defualt setting (lags = 1:12)")
        lags <- 1:12
    }
    else if (!all(base::round(lags) == lags)) {
        stop("Some of the inputs of the 'lags' argument are not integer type")
    }
    if (!is.numeric(margin)) {
        warning("The 'margin' parameter is not valid, using the defualt setting (margin = 0.2)")
        margin <- 0.2
    }
    if (!is.logical(Xshare)) {
        warning("The 'Xshare' parameter is not valid, please use only boolean operators.", 
            " Using the defualt setting setting (Xshare = TRUE")
        Xshare <- TRUE
    }
    if (!is.logical(Yshare)) {
        warning("The 'Yshare' parameter is not valid, please use only boolean operators.", 
            " Using the defualt setting setting (Yshare = TRUE")
        Yshare <- TRUE
    }
    if (stats::is.ts(ts.obj)) {
        if (stats::is.mts(ts.obj)) {
            warning("The 'ts.obj' has multiple columns, only the first column will be plot")
            ts.obj <- ts.obj[, 1]
        }
        df <- base::data.frame(time = stats::time(ts.obj), y = base::as.numeric(ts.obj)) %>% 
            dplyr::arrange(time)
    }
    else if (xts::is.xts(ts.obj) | zoo::is.zoo(ts.obj)) {
        if (!is.null(base::dim(ts.obj))) {
            if (base::dim(ts.obj)[2] > 1) {
                warning("The 'ts.obj' has multiple columns, only the first column will be plot")
                ts.obj <- ts.obj[, 1]
            }
        }
        df <- base::data.frame(time = zoo::index(ts.obj), y = base::as.numeric(ts.obj)) %>% 
            dplyr::arrange(time)
    }
    else {
        stop("The input object is not valid (must be 'ts', 'xts', or 'zoo')")
    }
    p_list <- lapply(base::seq_along(lags), function(i) {
        plotly::plot_ly(x = df$y %>% dplyr::lag(lags[i]), y = df$y, 
            type = "scatter", mode = "markers") %>% plotly::layout(xaxis = list(title = "", 
            range = c(base::min(stats::na.omit(df$y)), base::max(stats::na.omit(df$y)))), 
            yaxis = list(range = c(base::min(stats::na.omit(df$y)), 
                base::max(stats::na.omit(df$y)))), annotations = list(text = paste("Lag", 
                lags[i], sep = " "), xref = "paper", yref = "paper", 
                yanchor = "bottom", xanchor = "center", align = "center", 
                x = 0.5, y = 0.9, showarrow = FALSE))
    })
    p <- base::suppressWarnings(plotly::subplot(p_list, nrows = base::ceiling(base::length(p_list)/n_plots), 
        margin = margin, shareX = Xshare, shareY = Yshare) %>% 
        plotly::layout(title = main) %>% plotly::hide_legend())
    return(p)
} # In order to add personal title

# Load Data

Sales <- read.csv('/Users/lorenzoleoni/Desktop/Materiale Personale/Projects/Time Series Analysis - Car Sales In Italy/data/SalesITA.csv',
                  header = TRUE,
                  sep = ";",
                  dec = ",")

DATA PRE-PROCESSING

# Preparation of Data

t <- transpose(Sales)

names(t) <- Sales[,'Make'] # save columns names as models

t <- t[-c(1),] # drop row Make

t <- as.data.frame(lapply(t,as.integer)) # convert value into integers

row.names(t) <- names(Sales[,2:124]) # save rows names as periods

t$ITALY <- rowSums(t)

# Italy Sales Time Series

ITA <- ts(t$ITALY,
        start = c(2012, 1), 
        end = c(2022,03),
        freq = 12)

TIME SERIES ANALYSIS

# Summary Statistics

summary(ITA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4119  109093  133365  132127  155404  225448 

# Time Series

ts_plot(ITA,
     Xgrid = T, 
     Ygrid = T,
     Ytitle = "Units",
     Xtitle = "Year",
     title = "Time Series - Italy Sales",
     width = 2,
     color = "black")

NA

# Seasonal Italy Sales

ts_seasonal(ITA,
            type = "all",
            title = "Seasonal - Italy Sales")

NA

# Cycle and Seasonal Analysis

ts_heatmap(ITA,
           title = "Heatmap - Italy Sales")

NA

# Method of coefficient of variation of seasonal differences and quotients

i = 13
k = 0

D = vector() # Seasonal Difference
Q = vector() # Seasonal Quotient

while (k<length(ITA)-12){
  
  D[i-12] <- ITA[i]-ITA[i-12]
  Q[i-12] <- ITA[i]/ITA[i-12]
  
  i = i + 1
  k = k + 1
  
}

cv_d = sd(D)/mean(D) # coefficient of variation of the seasonal differences

cv_q = sd(Q)/mean(Q) # coefficient of variation of the seasonal quotients

dec = list()

seasonal_test <- mseastest(ITA, type = "pearson", m = 12)

if (abs(cv_q) <= abs(cv_d) & seasonal_test$is.multiplicative == T){
  print("The time series is multiplicative")
  dec = "mult"} else if (abs(cv_q) > abs(cv_d) & seasonal_test$is.multiplicative == F){
    print("The time series is additive") 
    dec = "add"} else{
      print("Tests have different results")
      }

[1] "The time series is multiplicative"

# Decomposition Time Series

ITA.dec <- decompose(ITA,
                   type = dec)

my_plot.decomposed.ts(ITA.dec, 
                      title = "Decomposition Time Series - Italy Sales",
                      cex.main = 1,
                      font.main = 1,
                      col.main = "#444444",
                      col.lab = "#444444",
                      cex.lab = 0.8,
                      cex.axis = 1,
                      xlab = "Year",
                      lwd = 2)

# Autocorrelation

ac <- acf(ITA,
          lag.max = 24,
          plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)
plot(ac, main = "")
title(main = "Autocorrelation Italy Sales")

# Partial Autocorrelation

pac <- pacf(ITA,
            lag.max = 24,
            plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)
plot(pac, main = "")
title(main = "Partial Autocorrelation Italy Sales")

# Seasonal Adjustment 

ITA_des <- diff(ITA,
              lag = 12)

ts_plot(ITA_des,
        Xgrid = T, 
        Ygrid = T,
        Xtitle = "Year",
        title = "TS Deseasonalized - Italy Sales",
        width = 2,
        color = "black")

NA

# Autocorrelation TS Deseasonalized

ac <- acf(ITA_des,
          lag.max = 24,
          plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)

plot(ac, main = "")
title(main = "TS Deseasonalized - Autocorrelation Italy Sales")

# Partial Autocorrelation TS Deseasonalized

ac <- pacf(ITA_des,
           lag.max = 24,
           plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)

plot(ac, main = "")
title(main = "TS Deseasonalized - Partial Autocorrelation Italy Sales")

# Trend Adjustment 

ITA_det <- diff(ITA,
              lag = 1)

ts_plot(ITA_det,
        Xgrid = T, 
        Ygrid = T,
        Xtitle = "Year",
        title = "TS Detrendalized - Italy Sales",
        width = 2,
        color = "black")

NA

# Autocorrelation TS Detrendalized

ac <- acf(ITA_det,
          lag.max = 24,
          plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)

plot(ac, main = "")
title(main = "TS Detrendalized - Autocorrelation Italy Sales")

# Partial Autocorrelation TS Detrendalized

ac <- pacf(ITA_det,
           lag.max = 24,
           plot = F)

par(cex.main = 1,
    font.main = 1,
    col.main = "#444444",
    font.lab=1,
    col.lab = "#444444",
    cex.lab = 0.8,
    cex.axis = 0.8)

plot(ac, main = "")
title(main = "TS Detrendalized - Partial Autocorrelation Italy Sales")

# Lag Plots

my_plot.lags.ts(ITA,
                main = "Lag Plots - Italy Sales")

NA

TIME SERIES MODEL SELECTION

HWES (Holt-Winter’s Exponential Smoothing Model)

# HWES - Auto

train.ts <- window(ITA,
                   start = c(2012,1),
                   end = c(2021,3))
test.ts <- window(ITA,
                  start = c(2021,3),
                  end = c(2022,3))

hw <- ets(train.ts,
          model = "ZZZ",
          restrict = F,
          allow.multiplicative.trend = T) # see plot seasonal in decompose

hw$method 

[1] "ETS(A,N,A)"

# Introductory Time Series With R - P. S. P. Cowperwait & A. V. Metcalfe (2009), pag 86

# Wavk Test for Trend Detection

notrend_test(train.ts,
             test = "WAVK",
             ar.method = "burg",
             ar.order = 2,
             BIC = F) 

    Sieve-bootstrap WAVK trend test

data:  train.ts
WAVK test statistic = 16.59, moving window = 11, p-value = 0.159
alternative hypothesis: (non-)monotonic trend.
sample estimates:
$AR_order
[1] 2

$AR_coefficients
      phi_1       phi_2 
0.545853992 0.000907434 

# Why Yule-Walker Should Not Be Used For Autoregressive Modelling.pdf
# funtimes.pdf, pag 25

# HWES - Additive Error

hw <- ets(train.ts,
          model = "ANM", # M seasonal how shown in previous section
          restrict = F)

hw.pred <- forecast(hw,
                    h = 12,
                    level = 0)

Train_Set <- train.ts

Test_Set <- test.ts

HW_Train <- hw$fitted

HW_Test <- ts(c(hw$fitted[111], hw.pred$mean), 
                    start = c(2021,3), 
                    end = c(2022,3), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Test_Set,HW_Train,HW_Test),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "HWES (A,N,M)",
        width = 2)

NA

# Metrics for HWES - Additive Error

print(paste("HWES Model with additive error --> AIC:",
            round(hw$aic,2),
            "BIC:",
            round(hw$bic,2),
            "AICc:",
            round(hw$aicc,2)))

[1] "HWES Model with additive error --> AIC: 2759.94 BIC: 2800.58 AICc: 2764.99"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(hw,
               lag = 24)

    Ljung-Box test

data:  Residuals from ETS(A,N,M)
Q* = 30.464, df = 10, p-value = 0.0007189

Model df: 14.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(hw$residuals)

    Shapiro-Wilk normality test

data:  hw$residuals
W = 0.71812, p-value = 2.663e-13

# HWES - Multiplicative Error

hw <- ets(train.ts,
          model = "MNM",
          restrict = F) # see plot seasonal in decompose

hw.pred <- forecast(hw,
                    h = 12,
                    level = 0)

HW_Train <- hw$fitted

HW_Test <- ts(c(hw$fitted[111], hw.pred$mean), 
                    start = c(2021,3), 
                    end = c(2022,3), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Test_Set,HW_Train,HW_Test),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "HWES (M,N,M)",
        width = 2)

NA

# Metrics for HWES - Multiplicative Error

print(paste("HWES Model with multiplicative error --> AIC:",
            round(hw$aic,2),
            "BIC:",
            round(hw$bic,2),
            "AICc:",
            round(hw$aicc,2)))

[1] "HWES Model with multiplicative error --> AIC: 2755.71 BIC: 2796.35 AICc: 2760.76"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(hw,
               lag = 24)

    Ljung-Box test

data:  Residuals from ETS(M,N,M)
Q* = 67.708, df = 10, p-value = 1.226e-10

Model df: 14.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(hw$residuals)

    Shapiro-Wilk normality test

data:  hw$residuals
W = 0.67974, p-value = 3.154e-14

SARIMA (Seasonal Autoregressive Integrated Moving Average)

# SARIMA Models

train.ts <- window(ITA, 
                   start = c(2012,1),
                   end = c(2021,3))
valid.ts <- window(ITA, 
                   start = c(2021,3), 
                   end = c(2021,9))
test.ts <- window(ITA, 
                  start = c(2021,9), 
                  end = c(2022,3))

sarima <- auto.arima(train.ts, 
                     seasonal = T) # (0,1,0)(2,0,0)

print(paste("Auto SARIMA Model --> AIC:",
            round(sarima$aic,2),
            "BIC:",
            round(sarima$bic,2)))

[1] "Auto SARIMA Model --> AIC: 2544.37 BIC: 2552.47"

sarima <- Arima(train.ts, 
                order = c(2,1,3), 
                seasonal = c(1,1,1))

print(paste("SARIMA Model based on theory --> AIC:",
            round(AIC(sarima),2),
            "BIC:",
            round(BIC(sarima),2)))

[1] "SARIMA Model based on theory --> AIC: 2252.2 BIC: 2272.88"

sarima <- Arima(train.ts, 
                order = c(0,1,0),
                seasonal = c(1,1,1))

print(paste("SARIMA Model near to auto.arima --> AIC:",
            round(AIC(sarima),2),
            "BIC:",
            round(BIC(sarima),2)))

[1] "SARIMA Model near to auto.arima --> AIC: 2256.41 BIC: 2264.16"

# Auto SARIMA Model

sarima <- auto.arima(train.ts, 
                     seasonal = T)

sarima.pred <- forecast(sarima,
                          h = 6,
                          level = 0)

Train_Set <- train.ts

Validation_Set <- valid.ts

Auto_SARIMA_Train <- sarima.pred$fitted

Auto_SARIMA_Validation <- ts(c(sarima.pred$fitted[111], sarima.pred$mean), 
                    start = c(2021,3), 
                    end = c(2021,9), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Validation_Set,Auto_SARIMA_Train,Auto_SARIMA_Validation),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "Auto SARIMA(0,1,0)x(2,0,0)",
        width = 2)

NA

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarima,
               lag = 24)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)(2,0,0)[12]
Q* = 22.011, df = 22, p-value = 0.4593

Model df: 2.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarima$residuals)

    Shapiro-Wilk normality test

data:  sarima$residuals
W = 0.79083, p-value = 2.855e-11

# Scoring Auto SARIMA Model

ac <- accuracy(sarima.pred, valid.ts)

row.names(ac)[2] <- "Validation Set"

ac

                        ME     RMSE      MAE        MPE     MAPE      MASE       ACF1
Training set      599.8114 23462.75 14247.40  -4.595007 15.48494 0.7895236 -0.0542865
Validation Set -65053.2068 72025.55 65053.21 -65.816179 65.81618 3.6049414  0.2716382
               Theil's U
Training set          NA
Validation Set  2.432375

# SARIMA Model based on Theory

sarima <- Arima(train.ts, 
                order = c(2,1,3), 
                seasonal = c(1,1,1))

sarima.pred <- forecast(sarima,
                        h = 6,
                        level = 0)

Train_Set <- train.ts

Validation_Set <- valid.ts

SARIMA_Theory_Train <- sarima.pred$fitted

SARIMA_Theory_Validation <- ts(c(sarima.pred$fitted[111], sarima.pred$mean), 
                    start = c(2021,3), 
                    end = c(2021,9), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Validation_Set,SARIMA_Theory_Train,SARIMA_Theory_Validation),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMA(2,1,3)x(1,1,1)",
        width = 2)

NA

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarima,
               lag = 24)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)(1,1,1)[12]
Q* = 10.566, df = 17, p-value = 0.8782

Model df: 7.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarima$residuals)

    Shapiro-Wilk normality test

data:  sarima$residuals
W = 0.59844, p-value = 5.915e-16

# Scoring SARIMA Model based on Theory

ac <- accuracy(sarima.pred, valid.ts)

row.names(ac)[2] <- "Validation Set"

ac

                        ME     RMSE       MAE       MPE     MAPE      MASE          ACF1
Training set      366.7765 18151.56  8995.061  -8.79325 15.26784 0.4984638 -0.0004354358
Validation Set -14139.2749 19587.28 17037.117 -13.45274 15.46243 0.9441165 -0.4827313486
               Theil's U
Training set          NA
Validation Set 0.7041307

# SARIMA Model near to Auto

sarima <- Arima(train.ts, 
                order = c(0,1,0),
                seasonal = c(1,1,1))

sarima.pred <- forecast(sarima,
                          h = 6,
                          level = 0)

Train_Set <- train.ts

Validation_Set <- valid.ts

SARIMA_Near_Train <- sarima.pred$fitted

SARIMA_Near_Validation <- ts(c(sarima.pred$fitted[111], sarima.pred$mean), 
                    start = c(2021,3), 
                    end = c(2021,9), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Validation_Set,SARIMA_Near_Train,SARIMA_Near_Validation),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMA(0,1,0)x(1,1,1)",
        width = 2)

NA

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarima,
               lag = 24)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)(1,1,1)[12]
Q* = 20.262, df = 22, p-value = 0.5667

Model df: 2.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarima$residuals)

    Shapiro-Wilk normality test

data:  sarima$residuals
W = 0.65812, p-value = 1.027e-14

# Scoring SARIMA Model Near to Auto

ac <- accuracy(sarima.pred, valid.ts)

row.names(ac)[2] <- "Validation Set"

ac

                        ME     RMSE       MAE        MPE     MAPE     MASE       ACF1
Training set      609.2847 19638.63  9789.521  -3.439208 11.12877 0.542489 0.01234894
Validation Set -42767.0066 48579.44 43417.129 -42.306566 42.75743 2.405972 0.02946240
               Theil's U
Training set          NA
Validation Set  1.696081

SARIMAX (Seasonal Autoregressive Integrated Moving Average Exogenous)

# Load Data and MIB Time Series

MIB <- read.csv('/Users/lorenzoleoni/Desktop/Materiale Personale/Projects/Time Series Analysis - Car Sales In Italy/data/MIB_12-22.csv',
                  header = TRUE,
                  sep = ";",
                  dec = ".")

MIB <- ts(MIB$CLOSE, start = c(2012,1), frequency = 12)

ts_plot(MIB,
     Xgrid = T, 
     Ygrid = T,
     Ytitle = "Value",
     Xtitle = "Year",
     title = "Time Series - MIB",
     width = 2,
     color = "black")

NA

# SARIMAX with MIB

train.MIB <- window(MIB, 
                   start = c(2012,1),
                   end = c(2021,3))
valid.MIB <- window(MIB, 
                   start = c(2021,3), 
                   end = c(2021,9))

sarimax <- Arima(train.ts, 
                order = c(2,1,3),
                seasonal = list(order = c(1,1,1), period = 12),
                xreg=train.MIB)

sarimax.pred <- forecast(sarimax,
                          h = 6,
                          level = 0,
                          xreg=valid.MIB)

Train_Set <- train.ts

Validation_Set <- valid.ts

SARIMAX_MIB_Train <- sarimax.pred$fitted

SARIMAX_MIB_Validation <- ts(c(sarimax.pred$fitted[111], sarimax.pred$mean), 
                    start = c(2021,3), 
                    end = c(2021,9), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Validation_Set,SARIMAX_MIB_Train,SARIMAX_MIB_Validation),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMAX(2,1,3)x(1,1,1) - MIB",
        width = 2)

NA

# Metrics SARIMAX Model with MIB

print(paste("SARIMAX Model with MIB as regressor --> AIC:",
            round(AIC(sarimax),2),
            "BIC:",
            round(BIC(sarimax),2)))

[1] "SARIMAX Model with MIB as regressor --> AIC: 2243.7 BIC: 2266.96"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarimax,
               lag = 24)

    Ljung-Box test

data:  Residuals from Regression with ARIMA(2,1,3)(1,1,1)[12] errors
Q* = 9.455, df = 16, p-value = 0.8935

Model df: 8.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarimax$residuals)

    Shapiro-Wilk normality test

data:  sarimax$residuals
W = 0.75017, p-value = 1.86e-12

# Scoring SARIMAX Model with MIB

ac <- accuracy(sarimax.pred, valid.ts)

row.names(ac)[2] <- "Validation Set"

ac

                       ME     RMSE      MAE        MPE     MAPE      MASE          ACF1
Training set     1117.204 17405.24 10165.97  -7.257166 15.35067 0.5633501  0.0003806399
Validation Set -38564.501 44762.82 40622.73 -38.013023 39.44043 2.2511200 -0.0042804871
               Theil's U
Training set          NA
Validation Set  1.577647

# Load Data and Unemployment Rate Time Series

Unr <- read.csv('/Users/lorenzoleoni/Desktop/Materiale Personale/Projects/Time Series Analysis - Car Sales In Italy/data/Unemployment_Rate_12-22.csv',
                  header = TRUE,
                  sep = ";",
                  dec = ".")

Unr <- ts(Unr$Tax, start = c(2012,1), frequency = 12)

ts_plot(Unr,
     Xgrid = T, 
     Ygrid = T,
     Ytitle = "Value",
     Xtitle = "Year",
     title = "Time Series - Unemployment Rate",
     width = 2,
     color = "black")

NA

# SARIMAX with Unemployment Rate

train.Unr <- window(Unr, 
                   start = c(2012,1),
                   end = c(2021,3))
valid.Unr <- window(Unr, 
                   start = c(2021,3), 
                   end = c(2021,9))

sarimax <- Arima(train.ts, 
                order = c(2,1,3),
                seasonal = list(order = c(1,1,1), period = 12),
                xreg=train.Unr)

sarimax.pred <- forecast(sarimax,
                          h = 6,
                          level = 0,
                          xreg=valid.Unr)

Train_Set <- train.ts

Validation_Set <- valid.ts

SARIMAX_Unr_Train <- sarimax.pred$fitted

SARIMAX_Unr_Validation <- ts(c(sarimax.pred$fitted[111], sarimax.pred$mean), 
                    start = c(2021,3), 
                    end = c(2021,9), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Validation_Set,SARIMAX_Unr_Train,SARIMAX_Unr_Validation),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMAX(2,1,3)x(1,1,1) - Unemployment Rate",
        width = 2)

NA

# Metrics SARIMAX Model with Unr

print(paste("SARIMAX Model with Unr as regressor --> AIC:",
            round(AIC(sarimax),2),
            "BIC:",
            round(BIC(sarimax),2)))

[1] "SARIMAX Model with Unr as regressor --> AIC: 2233.05 BIC: 2256.31"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarimax,
               lag = 24)

    Ljung-Box test

data:  Residuals from Regression with ARIMA(2,1,3)(1,1,1)[12] errors
Q* = 14.107, df = 16, p-value = 0.5908

Model df: 8.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarimax$residuals)

    Shapiro-Wilk normality test

data:  sarimax$residuals
W = 0.85498, p-value = 4.905e-09

# Scoring SARIMAX Model with Unemployment Rate

ac <- accuracy(sarimax.pred, valid.ts)

row.names(ac)[2] <- "Validation Set"

ac

                       ME     RMSE      MAE        MPE     MAPE      MASE        ACF1
Training set     2101.942 15945.15 10607.78  -4.461201 13.61187 0.5878331  0.01243276
Validation Set -15993.636 19667.10 18798.90 -16.373564 18.31904 1.0417463 -0.20415951
               Theil's U
Training set          NA
Validation Set 0.5908627

TIME SERIES FORECASTING

# Train and Test Sets

train.ts <- window(ITA, 
                   start = c(2012,1),
                   end = c(2021,9))

test.ts <- window(ITA, 
                  start = c(2021,9), 
                  end = c(2022,3))

# SARIMA Model based on Theory

sarima <- Arima(train.ts, 
                order = c(2,1,3), 
                seasonal = c(1,1,1))

sarima.pred <- forecast(sarima,
                        h = 6,
                        level = 0)

Train_Set <- train.ts

Test_Set <- test.ts

SARIMA_Theory_Train <- sarima.pred$fitted

SARIMA_Theory_Test <- ts(c(sarima.pred$fitted[117], sarima.pred$mean), 
                    start = c(2021,9), 
                    end = c(2022,3), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Test_Set,SARIMA_Theory_Train,SARIMA_Theory_Test),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMA(2,1,3)x(1,1,1)",
        width = 2)

NA

# Metrics for SARIMA Model based on theory

print(paste("SARIMA Model based on theory --> AIC:",
            round(AIC(sarima),2),
            "BIC:",
            round(BIC(sarima),2)))

[1] "SARIMA Model based on theory --> AIC: 2385.19 BIC: 2406.35"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarima,
               lag = 24)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)(1,1,1)[12]
Q* = 19.666, df = 17, p-value = 0.2916

Model df: 7.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarima$residuals)

    Shapiro-Wilk normality test

data:  sarima$residuals
W = 0.63702, p-value = 1.37e-15

# Scoring SARIMA Model based on Theory

accuracy(sarima.pred, test.ts)

                     ME      RMSE      MAE       MPE      MAPE      MASE        ACF1
Training set   439.9404 17628.486 9329.957 -9.279028 16.367974 0.4693342 0.004670013
Test set     -4302.7059  9882.253 8604.020 -3.498895  8.143556 0.4328166 0.475375284
             Theil's U
Training set        NA
Test set     0.6970851

# SARIMAX with Unemployment Rate

train.Unr <- window(Unr, 
                   start = c(2012,1),
                   end = c(2021,9))
test.Unr <- window(Unr, 
                   start = c(2021,9), 
                   end = c(2022,2))

sarimax <- Arima(train.ts, 
                order = c(2,1,3),
                seasonal = list(order = c(1,1,1), period = 12),
                xreg=train.Unr)

sarimax.pred <- forecast(sarimax,
                          h = 6,
                          level = 0,
                          xreg=test.Unr)

SARIMAX_Unr_Train <- sarimax.pred$fitted

SARIMAX_Unr_Test <- ts(c(sarimax.pred$fitted[117], sarimax.pred$mean), 
                    start = c(2021,9), 
                    end = c(2022,3), 
                    frequency = 12)

ts_plot(cbind(Train_Set,Test_Set,SARIMAX_Unr_Train,SARIMAX_Unr_Test),
        slider = T,
        Xgrid = T, 
        Ygrid = T,
        Ytitle = "Units",
        Xtitle = "Year",
        title = "SARIMAX(2,1,3)x(1,1,1) - Unemployment Rate",
        width = 2)

NA

# Metrics SARIMAX Model with Unr

print(paste("SARIMAX Model with Unr as regressor --> AIC:",
            round(AIC(sarimax),2),
            "BIC:",
            round(BIC(sarimax),2)))

[1] "SARIMAX Model with Unr as regressor --> AIC: 2368.93 BIC: 2392.73"

# Plot Residuals and Ljung-Box Test - H0: The level of correlation between the series and its lag is equal to zero, and therefore the series observations are independent

checkresiduals(sarimax,
               lag = 24)

    Ljung-Box test

data:  Residuals from Regression with ARIMA(2,1,3)(1,1,1)[12] errors
Q* = 16.857, df = 16, p-value = 0.3949

Model df: 8.   Total lags used: 24

# Shapiro-Wilk Test - H0: Normally Distributed

shapiro.test(sarimax$residuals)

    Shapiro-Wilk normality test

data:  sarimax$residuals
W = 0.77396, p-value = 3.891e-12

# Scoring SARIMAX with Unemployment Rate

accuracy(sarimax.pred, test.ts)

                    ME      RMSE      MAE       MPE      MAPE      MASE         ACF1
Training set  1920.983 16486.468 10250.30 -4.861421 13.522537 0.5156313 -0.015423834
Test set     -4809.261  6780.918  6551.21 -4.811418  6.447388 0.3295521  0.002396678
             Theil's U
Training set        NA
Test set     0.4666134

# Delete last variables

rm(seasonal_test,
   cv_d,
   cv_q,
   MIB,
   Test_Set,
   test.ts,
   Train_Set,
   train.ts,
   valid.ts,
   Validation_Set,
   ITA.dec,
   t,
   SARIMAX_Unr_Train,
   SARIMAX_Unr_Test,
   train.Unr,
   test.Unr,
   valid.Unr,
   SARIMA_Theory_Train,
   SARIMA_Theory_Validation,
   SARIMA_Theory_Test,
   SARIMAX_Unr_Validation,
   SARIMAX_MIB_Train,
   SARIMAX_MIB_Validation,
   train.MIB,
   valid.MIB,
   Auto_SARIMA_Train,
   Auto_SARIMA_Validation,
   SARIMA_Near_Train,
   SARIMA_Near_Validation,
   hw,
   hw.pred,
   HW_Train,
   HW_Test,
   ac,
   i,
   k,
   D,
   Q,
   pac,
   my_plot.decomposed.ts,
   my_plot.lags.ts,
   ITA_det,
   ITA_des)